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ABSTRACT 

We use the optical to mid-infrared coverage of the NEWFIRM Medium-Band Survey (NMBS) to char- 
acterize, for the first time, the properties of a mass-complete sample of 14 galaxies at 3.0 < z < 4.0 with 
M stal > 2.5 x 10 11 M Q , and to derive significantly more accurate measurements of the high-mass end of the 
stellar mass function (SMF) of galaxies at 3.0 < z < 4.0. The accurate photometric redshifts and well-sampled 
SEDs provided by the NMBS combined with the large surveyed area result in significantly reduced contribu- 
tions from photometric redshift errors and cosmic variance to the total error budget of the SMF. The typical 
very massive galaxy at 3.0 < z < 4.0 is red and faint in the observer's optical, with a median r-band magni- 
tude of (r tot ) =26.1, and median rest-frame U — V colors of (U-V) = 1.6. About 60% of the mass-complete 
sample have optical colors satisfying either the U- or the B-dropout color criteria, although ^50% of these 
galaxies have r > 25.5. We find that ^30% of the sample has SFRs from SED modeling consistent with zero, 
although SFRs of up to ~ 1 - 18 M Q yr" 1 are also allowed within 1 a. However, > 80% of the sample is 
detected at 24 /im, resulting in total infrared luminosities in the range (0.5-4.0) xlO 13 Lq. This implies the 
presence of either dust-enshrouded starburst activity (with SFRs of 600-4300 M© yr -1 ) and/or highly-obscured 
active galactic nuclei (AGN). The contribution of galaxies with M star > 2.5 x 10 11 M to the total stellar mass 
budget at 3.0 < z < 4.0 is ~ 8^ 3 %. Compared to recent estimates of the stellar mass density in galaxies with 
A^star « 10 9 - 10 11 Mq at z ~ 5 and z ~ 6, we find an evolution by a factor of 2-7 and 3-22 from z ~ 5 and z ~ 6, 
respectively, to z = 3.5. The previously found disagreement at the high-mass end between observed and model- 
predicted SMFs is now significant at the 3 a level when only random uncertainties are considered. However, 
systematic uncertainties dominate the total error budget, with errors up to a factor of ~ 8 in the densities at the 
high-mass end, bringing the observed SMF in marginal agreement with the predicted SMF. Additional system- 
atic uncertainties on the high-mass end could be potentially introduced by either 1) the intense star-formation 
and/or the very common AGN activities as inferred from the MIPS 24 fim detections, and/or 2) contamination 
by a significant population of massive, old, and dusty galaxies at z ~ 2.6. 

Subject headings: cosmology: observations — galaxies: evolution — galaxies: formation — galaxies: funda- 
mental parameters — galaxies: high-redshift — galaxies: luminosity function, mass func- 
tion — galaxies: stellar content — infrared: galaxies 



1. INTRODUCTION 

Understanding the formation mechanisms and evolution 
with cosmic time of galaxies is one of the major goals of ob- 
servational cosmology. An effective approach to understand 
the physical processes governing the assembly of galaxies 
(and their relative importance as a function of cosmic time) 
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is to directly measure the growth of the stellar mass in galax- 
ies. Galaxies can build their stellar mass both from in-situ star 
formation and/or merger events. The mean space density of 
galaxies per unit stellar mass, or stellar mass function (SMF), 
is one of the most fundamental of all cosmological observ- 
ables. The shape of the SMF retains the imprint of galaxy for- 
mation and evolution processes. Therefore, the SMF and its 
evolution with cosmic time represent a powerful tool to con- 
strain the physical mechanisms regulating the assembly and 
the evolution of galaxies. 

In the past decade, significant observational progress has 
been made in the measurement of the SMF of galaxies 
and its evolution with redshift. Using photometric redshifts 
derived from multi-waveband imaging surveys, measure- 
ments of the SMF o f galaxies are now rou t inely performed 
out to z ~ 5 (e.g . . iDickinson et alj 120031: IConselice et aTJ 
20051: iDrorv et all l2005t iFontana et all 120061: lElsner et al l 
20081 iPerez-Gonzalez et all 120081: iKaiisawa et al.1 120091: 
Marchesin i et al.l2009l) . The general consensus is that at z > 1 
the stellar mass assembly proceeds much more quickly than 
at lower redshifts. In particular, very recent measurements 
at z < 4 show a dramatic evolution of the SMF of galax- 
ies with redshift as well as evidence of mass-dependent evo- 
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lution of the SMF, with the low-mass end evolving more 
rapid l y than the high- mass end (i.e., iPerez-Gonzalez et al.l 
l2008tlMarchesini et alj|2009t) . 

Measurements of the SMF h ave been extended also to even 
larger redshifts (z ~4-7; e.g., iMcLure et alj|2009 l : iStark et alj 
2009), providing estimates of the stellar mass content of the 
universe when it was only ^800 Myr old. However, most of 
these studies only target Lyman break galaxies (LBGs), hence 
resulting in a potentially biased view of the universe against 
massive and evolved galaxies at such high redshifts. Whereas 
the discovery of a population of very massive and evolved 
galaxies at z > 5 has now been claimed by several groups (e.g. , 
Yan et alJl2006t iRodighiero et al.ll2007t IWiklind et~atll2008t 
Manci ni et al convincing evidence for the existence of 

galaxi es with M sta r > 3 x 10 11 M at z > 4 is still missing 
(e.g.. iDunlop eTal]|2007l) . 

Closely related to this issue is the very intriguing find- 
ing that the number density of the most massive galaxies 
(A^tar > 3 x 10 11 Mp) seems to ev olve very little from z ~ 4 
to z ~ 1.5 (Marches ini et al.l l2009). suggesting that the most 
massive galaxies in the universe were mostly in place already 
at z ~ 3 .5, and implying potentially severe disagreements with 
the predictions from the latest generations of semi-analytic 
models of galaxy formation. 

However, uncertainties on the observed SMF are still 
large, especially at the high-mass end and at high redshifts 
(Marchesin i et al.ll2009l) . At z < 3, the SMF error budget is 
now almost entirely dominated by systematic uncertainties 
caused by the different SED-modeling assumptions adopted 
to derive stellar masses (e.g., stellar population synthesis 
models or initial mass function (IMF)). Progress in reduc- 
ing the impact of these systematic uncertainties necessarily 
requires better calibrations of the stellar mass estimates (e.g., 
through measurements of the dynamical masses from stud- 
ies of their kinematics). At 3 < z < 4, instead, the contribu- 
tions of photometric redshift errors, small-number statistics, 
and sample variance (due to the relatively small probed vol- 
umes) are still significant, and dominate the total error budget 
at the high-mass end of the SMF. 

In this paper we take advantage of the high-quality 
data from the ultra-violet to the mid-infrared (MIR) avail- 
able through the NEWFIR M Medium-Band Survey (NMBS; 
van Dokkum et al. 2009b) to derive more accurate measure- 
ments of the high-mass end of the SMF of galaxies at 3.0 < 
z < 4.0 by significantly reducing the impact of random uncer- 
tainties and to characterize, for the first time, the observed and 
rest-frame properties of a mass-selected sample of galaxies at 
3.0 < z < 4.0. These results are made possible by the com- 
bination of accurate photometric redshifts and well-sampled 
spectral energy distributions (SEDs) of ^-selected galaxies at 
z > 1.5 delivered by the medium-band near-infrared (NIR) 
filters of the NMBS, as well as by its large surveyed area 
(~0.5 square degree). 

This paper is structured as follows. In § [2] we present the 
mass-selected (M stal - > 2.5 x 10 11 M Q ) sample used to mea- 
sure the SMF of galaxies at 3.0 < z < 4.0; in §[3] we present 
the observed and rest-frame properties of the galaxies in the 
mass-selected sample. The stellar mass function and densities 
of massive galaxies at 3.0 < z < 4.0 are presented in §0 while 
in §|5]the systematic effects caused by systematic uncertain- 
ties in the photometric redshift estimates are quantified. The 
results are summarized in §|6] We assume S1m = 0.3, il\ = 0.7, 
and Ho = 70 km s" 1 Mpc" 1 throughout the paper. All magni- 



tudes are on the AB system. 

2. SAMPLE SELECTION 

2.1. The NEWFIRM Medium Band Survey 

The sample is selected from the NMBS, a moder- 
ately wide, moderately d eep near-infrared imaging survey 
( Ivan Dokkum et al.ll2009bT) . The survey uses the NEWFIRM 
camera on the Kitt Peak 4m telescope. The camera images a 
28'x28' field with four arrays with a native pixel size of 0".4. 
We developed a custom NIR filter system for NEWFIRM, 
comprised of five medium bandw idth filters in the wavelengt h 
range 1 /im-1.7 fim. As shown in van Dokku m et al.l (12009b). 
these filters pinpoint the Balmer/4000 A breaks of galax- 
ies at 1.5 < z < 3.5, providing accurate photometric redshifts 
and improved measurements of the stellar population param- 
eters. The survey tar gets two 28'x28' fi elds: a subsection of 
the COSMOS field dScoville et alj|2007l) and a field contain- 
ing part of the AEGIS strip dDavis et al.ll2007|). Field posi- 
tions an d other information are given in Ivan Dokkum et"aT] 
(2009b). Both fields have excellent supporting data, in- 
cluding extremely deep optical ugriz data from the CFHT 
Legacy SurvexF] and deep Spitzer IRAC and MIPS imaging 
dBarmbv et al.ll2008t [Sandersll2007h . Spitzer-IRAC and MIPS 
photometry have been adde d following the procedure de- 
scribed in lWuvts et alj d2007l) . which uses a source-fitting al- 
gorithm developed by I. Labbe et al. (2010, in preparation) es- 
pecially suited for heavily confused images for which a higher 
resolution prior (in this case the /f-band image) is availableFI 
Reduced CF HT mosaics were kindly provided to us b y the 
CARS team dErben et alj|2009l: iHildebrandt et alj|2009l) . Ad- 
ditionally, in the COSMOS field, we include deep Subaru 
image s with the B{Vir + i + z + broad-band filters (Cap ak et al.1 
120071) . Subaru images with 12 intermediate-band filters from 
427 n m to 827 nm, and the CFHT X"s-band image (Ilbert et al. 
l2009h . In both the COSMOS and AEGIS fields, GALEX pho- 
tometry in the FUV (150 nm) and NUV (225 nm) passbands 
were added. The NMBS adds six filters: J u J 2 , /?, H u H 2 , 
and K. Filt ers characteristics of the fiv e medium band filters 
are given in Ivan Dokkum et al. (2009b). 

The NMBS is an NOAO Survey Program, with 45 nights al- 
located over three semesters (2008A, 2008B, 2009 A). An ad- 
ditional 30 nights were allocated through a Yale-NOAO time 
trade. The full details of the reduction, source detection, and 
generation of the photometric catalogs will be described in 
K. E. Whitaker et al. (in prep.). In the present analysis, we 
use a .^-selected catalog based on all the data obtained over 
the three semesters. The AEGIS catalog contains 17 filters 
and the COSMOS catalog contains 35 filters (FUV -8 /zm). 
The images were convolved to the same point-spread function 
(PSF) before performing aperture photometry, so as to limit 
any bandpass-dependent effects. Follo wing previous studies 
dLabbeet all 120031: lOuadriet aD 120071) . the photometry was 
performed with SExtractor in relatively small "color" aper- 
tures which optimize the S/N ratio. Total magnitudes in each 
band were determined from an aperture correction computed 
from the K-banA. The aperture correction is a combination 

10 |http : / /www, cfht .hawaii . edu/Science/CFHTLS/| 
11 The IRAC fluxes measured in this work have been com- 
pared with the publicly available IRAC photometry over COSMOS 
I http ://irsa . ipa c ■ caltech . edu/ data/ COSMOS/tables/scosmos/l 
Ilbert et al. 2009) and AEGIS I http : / /www . cf a ■ harvard ■ edu/irac/egs/| 
Barmby et al. 2008). The agreement is excellent, with systematic differences 
of ~2%. 
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of the ratio of the flux in SExtractor's AUTO aperture to the 
flux in the color aperture and a point-source-based correction 
for flux outside of the AUTO apertur e, thereby enabling us to 
calculate total magnitudes (see, e.g.. lLabbe et al.ll2.Q03b . The 
/T-band completeness limit of the NMBS catalog adopted in 
this work is K = 23.15 mag. Stars are flagged based on their 
observed U-Jl and Jl -K colors, where the stellar sequence 
cleanly separates from the bulk of galaxies in color space (see 
K. E. Whitaker et al. 2010, in prep., for more details). The 
total number of objects in the ^-selected sample is 52259, 
27520 of which are in the COSMOS field. 

2.2. Photometric Redshifts 

Phot ometric redshifts were determined w ith the EAZY 
code dBrammer. van Dokkum. & Coppi 2008), using the full 
FUV — 8 pm spectral energy distribution (SEDs) (FUV-K 
for objects in the ^50% of our AEGIS field that does not 
have Spitzer coverage) and z max = 6.0 (the maximum allowed 
redshift within EAZY). For this study we have used the pho- 
tometric redshift z vea k output by EAZYP1 Publicly available 
redshifts in the COSMOS and AEGIS fields indicate that the 
redshift errors are very small at <r z /(l +z) < 0.02 at z spec < 1 .5. 
Specifically, the photometric redshifts in COSMOS are in ex- 
cellent agreement with the spectroscopic redshifts made pub - 
licly available by the zCOSMOS survey dLilly et all 120071) . 
with a z /(l+z) = 0.008 for 1444 objects at z spec < 1.5. We 
also find excellent agreement between the photometric and 
spectroscopic redshifts for a larger sample of 23 1 3 objects at 
Zsper. < 1-5 in AEGIS from the DEEP2 survey dDavis et alJ 
120031) . with cr z /(l +z) = 0.017. Both fields have very few 
catastrophic failures, with only 3% > 5 a outliers. Although 
there are very few spectroscopic redshifts of optically-faint 
^-selected galaxies in these fields, we note that we found 
a similarly small scatter (cr z /(l + z) < 0.02) in a p ilot pro- 
gram targeting galaxies fr om the iRriek et al. (1200 8) near-IR 
spectroscopic sample (see van Dokkum et al l2009bl) . Spec- 
troscopic redshifts al so exist for 125 LB Gs at z ~ 3 within 
the AEGIS field from lSteidel et alJ (120031) . for which we find 
c z /(l +z) = 0.045, with 10% > 5 a outliers. From the formal 
EAZY errors listed in Table [T] we find typical er z /(l+z) = 
0.04, perfectly consistent with the scatter ct z /(1 +z) found for 
LBGs. We conclude that, in the regime of interest in this pa- 
per, the errors of the photometric redshifts are larger than at 
Zspec < 1-5, as they are dominated by random errors in the 
photometry. 

The observed spectral energy distributions (SEDs) with 
best-fit EAZY templates over-plotted are shown in Figure [T] 
together with the EAZY redshift probability distributions. As 
shown in Figure Q] the medium-band filters Hi and Hi allow 
us to identify the redshifted Balmer/4000A breaks within the 
H band, improving the accuracy of the photometric redshift 
estimates with respect to previous analysis with only broad- 
band photometry. 

Rest-frame colors were measured using the best-fit EAZY 
templates, as described inlBrammer et alJ (12009) and, in par- 
ticular, in IWhitaker et alJ d2010t) . Briefly, from the best-fit 
EAZY template, we comp uted the rest-frame U—V colors fol- 
lowing the method used by lWolf etafl (120031) in the COMBO- 

12 The default template set used in this work consists of seven templates: 
the six templates taken from the optimized template set of EAZY, but aug- 
mented with emission lines, and a template of a 12.5 Gyr old single stellar 
population. In §[5]we consider the case of an additional template, consisting 
of a dust-obscured (Ay = 3 mag), old (1 Gyr) population. 



17 survey. We used the Mafz Apellaniz| (120061) filter defini- 
tions and used the direct template fluxes to determine U — V. 
When using closely-spaced medium-band observed filters, the 
template fluxes are found to be more robust than interpolating 
between observed filters (see Brammer et al. 2010, in prep.). 

2.3. SED Modeling 

Stellar masses and other stellar population parameters 
were determined with FAST (Krieketal. 2009), fixing the 
redshift to the EAZY output (or the spectroscopic red- 
shift when available). For consistency with previously 
published SMF measurements and straightforward compar- 
iso ns, we assumed t he d efault SED-modeling assumptions 
of iMarchesini et alJ (120091). i.e., s tellar populat ion synthe- 
sis mo dels of iBruzual & Charlotl (12003b . the ICalzetti et all 
(2000) reddening law with Ay values rangi ng from to 4 in 
step of 0.2 mag, solar metallicity, pseudo-Kroupa (2001^1 
initial mass function (IMF), and three star formation histo- 
ries (SFHs): a single stellar population (SSP), a constant 
star formation history (CSF), and an exponentially declin- 
ing SFH with an e-folding timescale of 300 Myr (t3oo). In 
order to quantify the systematic uncertainties due to differ- 
ent SED-modeling assumptions on the derived stellar popu- 
lation properties (i.e., M star , age, star formation rate, and Ay) 
of the 3.0 < z < 4.0 sample, we ha ve also assumed the stel- 
lar population synthesis models of Maraston (2005) with a 
Kroupa| (12001b IMF and exponentially declining star forma- 
tion histories with values of the e-folding timescale ranging 
from 100 Myr to 10 G yr in step of 0.2 dex. We refer to 
Marchesin i et alJ (12009b for a detailed analysis of the system- 
atic uncertainties on the SMF measurements due to the differ- 
ent SED-modeling assumptions. Figure [3] shows the observed 
SEDs of the mass-selected sample at 3.0 < z < 4.0 together 
with the best-fit stellar population models from FAST for our 
two sets of SED-modeling assumptions. 

2.4. The 3.0 < z < 4.0 Mass-selected Sample 

We constructed a mass-selected sample of galaxies at 3.0 < 
z < 4.0 to study their observed and rest-frame properties, as 
well as to derive more accurate measurements of the high- 
mass end of the SMF of galaxies at 3.0 < z < 4.0. 

The redshift-dependent completeness limit in stellar mass 
h as been estimated followi ng the approach described in detail 
in lMarchesini et afl (12009b . which exploits the availability of 
samples with different depths. The completeness of a sample 
is estimated empiric ally from the available deeper samples, 
name ly, the FIRES (Labbe et al. 2003; Forster Schreib er et alJ 
I2006T) and the FIREWORKS (IWuyts et alJ 12008b catalogs. 
Briefly, to estimate the redshift-dependent stellar mass com- 
pleteness limit of the NMBS sample, we first selected galaxies 
belonging to the available deeper samples. Then, we scaled 
their fluxes and stellar masses to match the /T-band complete- 
ness limit of the NMBS sample. The upper envelope of points 
in the (M star . sca i e d-z) space, encompassing 95% of the points, 
represents the most massive galaxies at the considered flux 
limit (K= 23.15 for the NMBS catalog adopted in this work), 
and so provides a redshift-dependent stellar mass complete- 
ness li mit for the NMBS sample. We refer to Marchesi ni et alJ 
(2009) for a detailed description of this method. The result- 
ing completeness in mass of the NMBS catalog used in this 

SED modeling was performed using a Salpete^ (1953) IMF with lower 
and upper mass cutoffs of 0.1 Mq and 100 Mq, an d the de r ived s tellar 
masses and star formation rates were scaled to a pseudo- Kroupa 1 2001) IMF 
by dividing by a factor of 1.6. 
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FIG. 1. — Observed SEDs of the mass-selected sample at 3.0 < z < 4.0. Filled circles are the observed fluxes, in units of 10~" erg cnr 2 s A -1 , with 
corresponding I <j errors. The blue symbols are the photometric points from NMBS. The solid gray curves represent the best-fit EAZY templates. The dark gray 
filled regions represent the EAZY redshift probability functions. The vertical red line is the adopted redshift from EAZY (z pea k, as specified in S I2.2K while the 
shaded gray regions are the 1 a allowed values for the photometric redshifts. Also listed are the NMBS identification number, the adopted photometric redshift 
from EAZY, and the observed total K- and r-band magnitudes. 
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TABLE 1 

Mass-selected sample of 3.0 < z < 4.0 galaxies 



ID 


r H 


K z 


logM star 


SFR 


A V 


log Age 




(mag) (mag) 


(mag) 


(M ) 


(M Q yr l ) 


(mag) 


(yr) 



CI -4890 

Cl-6110 

CI -7340 

Cl-15182 

Cl-15367 

Cl-18825 

Cl-19536 

Cl-21316 

CI -22857 

Cl-23152 

A2-6835 
A2-15753 
A2- 18070 
A2-24511 



26.27 

> 28.4 
28.87 
25.88 
28.49 
27.18 
25.49 
25.60 

> 28.3 
23.22 

26.19 
25.12 
25.35 
25.32 



24.90 
24.91 
24.32 
22.97 
23.86 
23.63 
22.66 
23.79 
24.64 
20.96 

23.61 
22.79 
23.04 
22.97 



22.39±0.16 
23.11±0.16 
23.07±0.15 
21.62±0.09 
23.07±0.19 
22.62±0.12 
21.65±0.06 
22.29±0.16 
23.09±0.19 
20.31±0.02 

22.28±0.13 
22.25±0.06 
22.33±0.11 
21.89±0.09 



3 47+0.34 
J - H '-0.50 

3.58J; 
3 41 +0.38 

-0.27 

3 56+*- 11 
J - ,D -0.I1 
3 73+0-22 
-'■'-'-0.06 

3 05+°' 19 

3 ,Q+0.0? 

3.54+fl.20 

°- Ay -Q.Q6 

3 07+ - 54 
J - u/ -0.70 
3 14+O.IO 
J - 1H -0.09 

3.08«; 
3.76tH 



11 70+°- 17 
H.41+0.07 

1 1 46+1i:07' 

1 1 54+004 
11 - ;IH -0.05 
11 71+005 

1 1 40+004 

I , <r 5 +0.03 

11.52; 001 

II 42+^ : ^ 

11.48** 

11.44« 
11.68+°;' J 



67 +o.07-> 



3q+o:o8 
a 1 +8:o? 

J1 -0.11 

45+o.os 



-0.09^ 
51 +0.12 



-0.02 



) 

2 6 +o:io-, 

-0.22 - 1 
99+0.09 

-0.03 > 
39+0. 11-) 



37 



48+ 



" 7 -0.06' 
30+0. 12x 



40.7+1^3 (33.9+J" ) 



8.3+' 1 / (2.3%) 
O.oJ (0.0+^) 
14.8+^-3 (13.2+1 8 /) 
1.9^-2 (2.0+$) 



316«f (251+$') 
0.0^(0.0^) 
0.1«;f(0.9«|) 

178^5 (H2+M5) 

148^| (60.3+™) 
166+™(151+ 7 ") 
170+™ (97.7+^) 



i 4+O.6 |"i 9+1 .2-, 
1.0+0-6 (0 8 +0.4 } 

l-2l;l(0.7 

1 



0.4 ' 

6 +o.2 , 



02 ' 
+0.4 , 



u -0.4 

i.4i 

0.8^1 
i-ol 

0.8* 
0.811 

2.(rt>j(1.7+°;§) 
1.4^;|(1.0+g) 
1.6l : |(1.5l|) 



d-i-Sij) 

(0.8l|) 

(o-i1:f) 

(0.4+° :l ,) 
(0.711) 



9j+o.i 
9 2+^ : ^ 
9 2+W 

o -j+0.2 

s - 3 -o.i 
9.2+ -0 

9 3+0:0 

1+0.1 
'-0.1 

9 i+o.i 

-14 
9.2+0-i 

8.1+ 



9.1+ 



(9.2+01) 
(9.2« : «) 
(9.2lj) 
(8.3+g : |) 
(9.l|T) 
(9.21-1) 

(9.o;o ; l) 

(9.1^') 

: (9.21; ) 
(8.0+0 J) 



9.3+0-0 (9 3 
9.3+o:o (8.8+0-=) 
q -i+O.O ,-q 9+0. l\ 
'•-0.1 ^* -Q.4' 
8.9+0;3 (8.4+0- 4 ) 



+0.0 
1.7' 
5-, 



NOTE. — "CI" and "A2" refer to the COSMOS and AEGI S fields respect ively. The listed redshift i s the a dopted best-fit EAZY redshift z peak . 
The st ellar population par ameters were derived u sing a pseudo- Kroupa 12001) IMF, Bruzual & Chariot (2003) stellar popul ation synthesis mo dels, 
and afCalzetti et al. 12000) extinction law (see j} |2.3> . Quoted errors are the 1 a confidence intervals output by FAST (see Kriek et al. 2009] for a 
detailed description of the adopted method in FAST to estimate c onfidence intervals). The values in parenthesis correspond to the best-fit stellar 
population parameters assuming a Kroupa 1 200 1) IMF, Maraston 12005) stellar population synthesis models, exponentially declining SFHs, and 
a Calzetti et al. (2000) extinction law (see § 12.3) . //-band magnitudes are derived by averaging the HI- and Z/2-band magnitudes. 
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FIG. 3. — Observed SEDs of the mass-selected sample of galaxies at 
3.0 < z < 4.0. Filled circles are the observed fluxes in arbitrary units, with 
corresponding 1 <r errors. The blue symbols are the photometric points from 
NMBS. The solid green curves represent t he best-fit FAST te mplates using 
the stellar population synthesis models of Bruzual & Chariot (2003). The 
solid red curves represent the best-fit FAST templates using the stellar popu- 
lation synthesis models of Maraston 1 2005). Both stellar population synthesis 
models provide generally good fits to the observed SEDs. 



work is M star = 10 1L4 ° M Q « 2.5 x 10 11 M over the targeted 
redshift range 3.0 < z < 4.0. 

The resulting mass-selected sample of galaxies at 3.0 < z < 
4.0 contains 14 sources with M star > 10 1L40 M Q (10 from the 
COSMOS field and 4 from the AEGIS field) over an effective 
area of 0.44 square degrees. The sample is listed in Table Q] 
along with the observed r-, H-, and /T-band total magnitudes, 
adopted EAZY best-fit redshifts and 1 a errors, and FAST 
best-fit Mstar. SFR, Ay, and age with corresponding 1 a errors. 

As shown in Table Q] the typical random error on the esti- 
mated stellar masses of the mass-selected sample is ~ 0.1 dex 
for the default set of SED-mod eling assumptions (which 
uses the iBruzual & Chariot! 120031 models), and ~ 0.16 dex 
for the other set (which adopts the [Maraston 2005 models). 
These errors are in good agreement with the errors on stellar 
mass due to photo metric redshift uncertainties estimated by 
iTavlor etail (120091) . with a typical error on the stellar mass of 
~ 0.1 dex for photometri c redshift er r ors of <x z /(l +z) = 0.035 
at z < 1.5. As shown by Taylo r et alj ( 120091) . in a photometric 
redshift survey, the stellar mass estimates are relatively robust 
to random photometric redshift errors, due to the similar (but 
opposite) systematic effects on luminosities and stellar mass- 
to-light ratios caused by random photometric redshift errors. 

Figure [4] shows the images of the mass-selected sample at 
3.0 < z < 4.0 in the different filters, from the M-band to the 
24 fim. Spitzer-MWSi channel. 

In order to exclude contamination of the mass-selected sam- 
ple due to blending, we have also inspected the higher-spatial 
resolution images available over COSMOS and AEGIS. For 
COSMOS, we have used the CFHT-WIRCAM K & -bwA im- 
age (FWHM-0.8") and the HST-ACS 7 F8 i4w-band imagesQ 
For AEGIS, we have used the HST-ACS /pg^wband im- 
agesQ In the HST images, only Cl-19536, Cl-23152, and 
A2-15753 are detected, whereas the other sources do not show 



14 Available at http : //irsa . ipac . caltech . edu/data/COSMOS/datasets . 

15 Available athttp://aegis.ucolick. org/acs_data_descrip . html 
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FIG. 4. — Images of the mass-selected sample of galaxies at 3.0 < z < 4.0. From left to right, the columns show the CFHTLS «, g, r, i, and z. images, the NMBS 
71, J2, 73, HI, H2, and K images, the Spitzer-IRAC 3.6 fim, 4.5 fim, 5.8 (im, and 8.0 pm images, and the Spitzer-MIPS 24 pm image. Each cutout is 6"x6" 
on a side. 



any obvious detection. All three detected sources appear to be 
resolved in the ACS images, indicating that the /-band fluxes 
are not dominated by a point-like component. Specifically, 
Cl-19536 and A2-15753 are extended and quite elongated. 
All the sources but C 1-23 152 appear very isolated in the ACS 
images, consistent with the ground-based images. The HST- 
ACS /F8i4w-band image of C 1-23 152 reveals two fainter knots 
at a distance of ~ 1.1". The two knots contribute about 11% 
to its total flux in the ACS image. Inspection of the As-band 
image over COSMOS reveals that all selected massive galax- 
ies are single, isolated objects, including Cl-23152, showing 
no obvious signature of the two knots. If the photometry of 
Cl-23152 is equally affected by the two knots at all wave- 
lengths, the shape of its SED would not be affected, and the 
resulting stellar mass would be smaller by ~ 0.05 dex, not 
changing the results of this paper. On the contrary, if the 
contribution of the two knots changes as a function of wave- 
length, the observed SED would change accordingly, mak- 
ing it harder to predict how the derived stellar mass would 
be affected. A rough estimate of this effect was derived by 



re-fitting the observed SED of Cl-23152 after assuming that 
only the optical fluxes are affected by the two knots. The 
resulting stellar mass is only ~ 0.03 dex smaller than that 
estimated with the current photometry, implying that the de- 
rived stellar mass for C 1 -23 1 52 is not likely to be significantly 
biased by the two knots. We therefore conclude that none 
of the observed objects seem to be affected by blending is- 
sues, which might have resulted in systematically biased stel- 
lar mass estimates. Higher spatial-resolution NIR imaging is 
however required to confirm this. 

Finally, we note that no a priori exclusion of active galac- 
tic nuclei (AGNs) has been performed in our mass-selected 
sample. 

3. PROPERTIES OF VERY MASSIVE GALAXIES AT 3.0 < Z < 4.0 

We use our mass-selected sample of 14 galaxies to de- 
termine the median and dispersion in observed and rest- 
frame properties of the most massive galaxies (M sta r > 2.5 x 
10 11 M ) at 3.0 < z < 4.0. Table |f] lists the median and 
25th/75th percentiles of the distributions of observed r-band 
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TABLE 2 
Observed and rest-frame 
properties of the 3.0 < z < 4.0 
mass-selected sample 



Quantity 


25% 


Median 


75% 


r tot (obs) 


25.3 


26.1 


27.8 


H-K (obs) 


0.75 


1.16 


1.42 


Vtot (rest) 


-24.2 


-23.5 


-23.3 


U-V (rest) 


1.26 


1.64 


1.90 


/3 (rest) 


-0.56 


-0.36 


0.07 



magnitude and H-K color; the rest-frame V-band magnitude 
and U-V color; and rest-frame UV slopes, parametrized by 
F\ oc A* 3 . Rest-frame UV slope s were determined from the 
best-fitting SEDs, following the lCalzetti et all (119941) method 
of fitting to the 10 rest-frame UV bins defined by those au- 
thors. 

Figure[5]shows the distributions of the observed H-K color 
(top panel), rest-frame U-V color (middle panel), and rest- 
frame UV slopes of the mass-selected sample (bottom panel). 
For comparison, we have also plotted 1) the distribution of 
rest-frame U-V colors and UV slopes of a mass-selected 
sampl e of galaxies at 2 < z < 3 with M star > 6 x 10 10 M Q 
from Ivan Dokkum et all (120061) (orange histogram); 2) the 
distribution of rest-frame U-V colors and UV slopes of 
the ga laxies that would be selec ted as LBGs from the sam- 
ple of van Do kkum et al.l d2006l) (purple histogram); and 3) 
the distribution of rest -frame UV slopes of a z ~ 3.7 sam- 
ple of galaxies from Bram mer & van Dokkuml (120071) se- 
lected with the color criterion H-K > 0.9 to have prominent 
Balmer/4000 A breaks between the H and K bands (cyan his- 
togram). 

The typical very massive galaxy at 3.0 < z < 4.0 (median 
stellar mass (M star ) ~3x 10 11 M ) is red and faint in the 
observer's optical, with (r tot ) = 26.1. Most galaxies (10 out 
of 14) would be selected as H- K galaxies with H-K > 0.9 
( Bram mer & van Dok kum 2007). Of the four galaxies with 
H-K < 0.9, only two galaxies have H-K color significantly 
smaller than the H-K criterion. This highlights the effi- 
ciency of this color technique in selecting galaxies at z > 3 
with prominent breaks in the rest-frame optical, although the 
fraction of interlopers selected by this color technique remains 
uncertain. 

From Table Q] 40%-60% of the very massive galaxies are 
characterized by ages consistent with the age of the universe 
at the targeted redshifts (~ 1.6-2.1 x 10 9 yr). About 30% of 
the very massive galaxies, namely Cl-6110, Cl-15182, Cl- 
22857, and Cl-23152 have SFR estimates from SED model- 
ing consistent with no star formation activity to within 1 a, in- 
dependent of the specific SED-modeling assumptions adopted 
in FAST. Of the remaining galaxies, 4 have very large SFRs, 
of the order of a few hundreds solar masses per year. We 
stress that the estimated ages and star formation rates from 
SED modeling are quite uncer tain, even with the h igh-quality 
dataset used in this work (e.g., Muzzin et al. 20091). 

3.1. Rest-frame UV 

The rest-frame U-V colors range from U-V = 1.01, typi- 
cal of nearby irregular gal axies, to U-V = 2.2, typical of lo- 
cal elliptical galaxies (e.g., Fukugi ta et alJI 1995b . The median 
(U-V) = 1 .64 is similar to local Sb spiral galaxies. As shown 
in the middle panel of Figure [3] the mass-selected sample of 
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FIG. 5. — Top panel: Distribution of observed H — K colors of the mass- 
selected sample of galaxies at 3.0 < z < 4.0 (solid black line); the blue 
hatched area represents the distribution for those galaxies that would be se- 
lected as either U- or fi-dropout galaxies based on their observed optical 
colors; the red filled area represents the distribution for those galaxies that 
would not be selected as either [/-or B-dropout galaxies based on their ob- 
served optica l colors; the vertical gr ay line represents the H — K criterion 
adopted in Brammer & van Dokkum 1 2007) to select galaxies at z ~ 3.7. 
Middle panel: Distribution of rest-frame U — V colors of the mass-selected 
sample; colors as in the top panel; the orange solid line represents the distri- 
butio n of rest-frame U - V colo rs of the mass-selected sample at 2 < z < 3 
from Ivan Dokkum et al. (2006); the purple solid line represents the distri- 
bution of rest-frame U — V colors of the galaxies in the mass-selected sam- 
ple at 2 < z < 3 from [van Dokkum et al. 1 2006) that would be selected as 
LBGs. Bottom panel: Distribution of rest-frame UV slopes of the mass- 
selected sample; colors as in the middle panel; the solid cyan line represents 
the distribution of UV slopes o f the H — K selected sample at z ~ 3.7 of 
Brammer & van Dokkum (2007). 



Ivan Do kkum et al. (2006]) at z = 2.5, which is complete in stel- 
lar mass down to ~ 6 x 10 10 M Q (a factor of ~ 5 less than 
our sample), encompasses the range in U-V colors spanned 
by our z = 3.5 sample, with a median U-V color bluer by 
~ 0.1 mag with respect to our mass-selected sample. 

The median UV slope f3 is (f3) = -0.36, indicating a rela- 
tively fiat spectrum in Fx- The distribution of j3, ranging from 
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Fig. 6.— Color selection of LBGs at z~3.0±0.3 ISteidel et alJ2003D and , 
in the insert, of B-dropout galaxies at z ~ 3.8 ± 0.3 iBouwens et alj |2009l) . 
Objects falling in the gray shaded regions would be selected as LBGs or B- 
dropout galaxies. Out of the 14 very massive galaxies at 3.0 < z. < 4.0, 3 
would be selected as LBGs and 5 as B-dropout galaxies based on their colors. 
However, only 4 out of 8 of the dropout galaxies have r < 25.5. 

= -0.95 to = 1.10, is broadly consiste nt with the distribu- 
tion of massive galaxies at 2 < z < 3 from lvan Dokkum et alj 
(2006). As shown in the bottom panel of Figure [5] the distri- 
bution of is instead very different from the distribution seen 
for H - K > 0.9 galaxies at z ~ 3.7, wh ich shows a peak at 
~ -2 (Brammer & van Dokkum 2007). The observed dis- 
tribution of for the mass-selected sample at 3.0 < z < 4.0 
is also very different from the distributions seen for UV- 
selected galaxies at z ~ 2.5 and z ~ 4, which peak at ~ -1 .6 
and ~ ~l-8, respectively (e.g., [A delberger & Steideli pOOOt 
Ouch i et alJl2004[lHamietal]l2008t IBouwens et al.ll2009l) . 

The intrinsically different rest-frame UV properties of the 
mass-selected sample at 3.0 < z < 4.0 studied in this work 
and the typical UV-selected galaxies at these redshifts (i.e., 
U- and B-dropout galaxies) is also clear from Figure|6l which 
shows the location of the massive galaxies at 3.0 < z < 4.0 
in the U n GR and B435V6O6Z850 diagra ms commonly used to se- 
lect [/-dropout galax ies (i.e.. LBGs: [Steidel et alJl2003h and 
B-dro pout galaxies dGiavalisco et al.l 12004; Bouwe ns et alj 
2009), respectively. The colors plotted in Figure [6] are syn- 
thetic colors integrated from the best-fit FAST templates. 
About 57% of the galaxies in the 3.0 < z < 4.0 mass-selected 
sample have colors that satisfy either the U- or the B-dropout 
color criteria (gray shaded area in Figure [6]). Of these, three 
would be selected as {/-dropouts, and five as B-dropouts, 
based on their observed optical colors. However, ^50% of 
these UV-selected galaxies are fainter than r tot = 25.5, which 
is the observed optical limit of typical spectroscopic samples 
of LBGs. While the r tot = 25.5 cut is not relevant to the inclu- 
sion of our galaxies in the photometric window, it is relevant 
when considering our objects in the context of pre-existing 
LBG surveys. 

The rest-frame SEDs of the mass-selected sample at 3.0 < 
Z < 4.0 are shown in Figure [7j together with the median rest- 
frame SED from the data (solid blue curve) and the median 
best-fit templates from FAST (green and red solid curves). 
Figure UJ clearly shows the strongly-suppressed emission and 
the flatness of the spectrum in Fx in the rest-frame UV, as 
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FIG. 7. — Rest-frame SEDs of the mass-selected sample at 3.0 < z < 4.0 
(yellow filled circles). The SEDs are normalized to the flux at A rcs t = 4000 A. 
The solid blue curve represents the running med ian of the 15 neighboring 
points. The best-fit FAST templates adopting the Bruzual & Chariot ( 2003) 
models are shown as dark gray solid curves. The green a nd red solid curves 
repres ent the me dian best-fi t FAST templates adopting the Bruzual & Charlotl 
120031) and the Maraston 12005) models, respectively. The cyan solid 
curve represents the median SE D of the H — K galaxies at z ~ 3.7 from 
IBrammer & van Dokkuml UMfo . 

well as the strong Balmer/4000 A breaks in the rest-frame 
optical for the typical very massive galaxy at z ~ 3.5. Also 
plotted is the median rest-frame SED of the H — A"-se lected 
sample at z ~ 3.7 from lBrammer & van Dokkuml (T2007). Fig- 
ure UJ clearly shows that the rest-frame optical SED of the 
H — ^-selected sample and our mass-selected sample are very 
similar, characterized by strong rest-frame optical breaks. In 
contrast, their rest-frame UV SEDs are very different. The 
H - ^-selected galaxies are characterized by very blue rest- 
frame UV-optical colors. On the contrary, our mass-selected 
galaxies are generally red also in the rest-frame UV. 

The significant differences in observed and rest-frame prop- 
erties between the H - /T-selec ted galaxies at z ~ 3.7 from 
(IBrammer & van Dok kum 2007) and the very massive galax- 
ies at 3.0 < z < 4.0 selected in our study are very interesting, 
as most of our galaxies would also be selected as H-K galax- 
ies. The simplest explanation for the observed differences is 
the very different regime in stellar masses probed by the two 
sampl es. The H-K galaxies from Bra mmer & van Dokkuml 
(2007) were selected from the FIRES survey over an area of 
~ 31 arcmin 2 , with a median stellar mass a factor of ~ 15 
smaller than the median stellar mass of our mass-selected 
sample. The lack of very massive galaxies in FIRES is sim- 
ply caused by its very small surveyed area, as an area of > 
100 arcmin 2 is required to find one single object as massive as 
our mass-selected galaxies. The different stellar mass regime 
probed by the two samples suggests that lower-mass galaxies 
are characterized by much bluer rest-frame UV-optical colors 
than the most massive galaxies at these redshifts. This is also 
supported by the different rest-frame UV properties between 
our mass-selected sample and typical UV-selected B- and 
V-dropout samples, which have stellar masses in the range 
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10 9 - 10 11 M (e.g. IShaolev eFaIll200Tt iMaedis et al.ll20loh . 
From the SED modeling analysis, it appears that the differ- 
ences in the rest-frame UV properties between the H-K- 
selected galaxies and our mass-selected sample could be pri- 
marily due to the different amount of dust in the two samples, 
as their typical stellar ages (~ 1-1.5 Gyr) are broadly con- 
sistent within the uncer tainties. More specifically, the H -K- 
selected galaxies from Bram mer & van Dok kum (2007|) are 
characterized by a median dust extinction of (Ay) ~ 0.2, a 
factor of about ~ 7 smaller than the median extinction of 
our mass-selected sample ((Ay) = 1.4 mag). Indeed, signif- 
icant amount of dust seems to be quite ubiquitous in massive 
galaxies at 3.0 < z < 4.0. The amount of dust in our mass- 
selected sample is significantly higher than the amount of dust 
in local massive galaxies and in massive galaxi es at z ~ 1-2 
(A v ss 0.2-0.3 mag, e.g., IWhitaker et al.ll2010l) , further con- 
tributing to suppressing the rest-frame UV light in these mas- 
sive galaxiesFl 

3.2. Spitzer-MIPS 24 fjm data 

We have used the publicly available observations at 24 /im 
from Spitzer-MIPS provided by S-COSMOS and FIDEL0 to 
further constrain the activity in the most massive galaxies at 
3.0 < z < 4.0. The measured 24 /im fluxes with the corre- 
sponding 1 a errors are listed in Table [3] For two galax- 
ies over AEGIS, namely A2-6835 and A2-24511, no MIPS- 
24 /im data are available. The 24 /im flux of CI- 18825 is 
not reliable due to blending issues from a nearby very bright 
24 /im source. The MIPS cutouts are shown in Figures H] 

In the sample of 1 1 massive galaxies at 3.0 < z < 4.0 with 
MIPS coverage and no blending issues, ^80% have a MIPS 
24 /im fluxes significant at > 3 a. This is broadly consis- 
tent with the high fraction of MIPS -detected sources in the 
sample of IRAC-selected massive (M star ~ 10 10 — 10 1 1 M^ ) 
galaxies at z > 3.5 over GOODS-North (Ma ncini et aLll2009h . 
The fraction of MIPS detected massive galaxies in our sample 
increases up to ^90% for a > 1 a detection. Only CI -22857 
is undetected at 24 /im. 

For the redshift range targeted here, the 24 /im band probes 
the rest-frame wavelengths from ~ 4.8/im to ~ 7.1/mi, which 
includes the 5.27/mi, 5.7/im, and 6.2/im emissi on features 
from poly cyclic aromatic hydrocarbons (PAHs) dDraine & Lil 
2007). MIR emission at these wavelengths could arise 
from warm/hot dust and PAH molecules, heated by either 
dust-enshrouded star formation or AGN. The 24 /im emis- 
sion is widel y used to estimate the SFRs in high-re dshift 
galaxies (e.g.. iRigbv et al]|2008t iPapovich et al.i r2007). The 
dust-enshrouded SFRs can be estimated by transforming to- 
tal infrared lumi nosities (Lir = L(8- 1000 pm)) into SFRs 
(lKennicuttlll998l) . 

To convert the 24 /im emission to a total IR luminosit y 
we followed the approach presented in IWuvts et al.l (2008). 
Specifically, we used the infrared SED s of star-forming galax- 
ies provided by iDale & Heloul (120021) . The template set al- 
lows us to quantify the IR/MIR flux ratio for different heat- 
ing levels of the interstellar environment, parameterized by 

16 We note that performing the SED-modeling forcing Ay = results in 
stellar masses smaller by — 0.2 dex. However, the resulting \ 2 of the best- 
fit models are significantly worse than the \- values corresponding to the 
modeling allowing for dust. Moreover, the MIPS detections strongly suggest 
the presence of significant amount of dust. Therefore, we conclude that the 
dust-free assumption is not a realistic assumption. 

17 |http : / /irsa ■ ipac ■ caltech . edu/ data/ SPITZER/FIDEL/I 



TABLE 3 

Spitzer-MIPS 24 flM FLUXES AND DERIVED PROPERTIES OF 
THE 3.0 < Z < 4.0 MASS -SELECTED SAMPLE 





ID 







i-IR 


SFR 






no 13 1 ^1 

Liu 1^01 


HU^ vr" 1 1 
L1V10 yi J 


CI -4890 


206.2±32.8 


3.3±0.5£jl) 


3611±574(;2™) 


Cl-6110 


116.4±28.0 


2.3±0.5c|i) 


2458±591(!!'«) 


CI -7340 


133.0±28.9 


1.9±0.4d|) 


2078±451(+JgJ2) 


Cl-15182 


78.6±25.7 


1.5±0.5d|) 


1614±5280 
4275±787(«f 6 ) 


Cl-15367 


167.3±30.8 


3.9±0.7(«|) 


Cl-18825" 


blended 


Cl-19536 


61.6±24.6 


0.5±0.2Cj") 
4.0±0.7Cf°) 


584±233( + |°°) 
4330±765C}gg) 


Cl-21316 


177.1±31.3 


CI -22857 


<19.9 


<0.4 


<395 


Cl-23152 


110.5±27.6 


1.2±0.3(^) 


1331±333(^) 


A2-6835* 








A2-15753 


165.7±22.8 


1.3±0.3C5j) 
0.9±0.3O 


1371±378(+|57 } 


A2- 18070 


127.0±21.3 


918±308O 


A2-24511* 









NOTE. — "no reliable MIPS 24/im flux could be obtained for 
Cl-18825 due to blending issues; *no MIPS 24/^m data available 
for A2-6835 and A2-23152. The errors listed for L lR and SFR are 
computed using just the 24 fim photometric errors (values not in 
parenthesis) and the combination of the 24 fim photometric errors 
and the photometric redshift errors (values in parenthesis). 

dM(U) ~ U~ a dU, where M(U) represents the dust mass 
heated by an intensity U of the interstellar field. We com- 
puted the total i nfrare d luminosity Lir Q for each object for all 
IDale & Heloul (|2002) templates within the reasonable range 
from a = 1 for active galaxies to a = 2.5 for quiescent galax- 
ies. The mean of the resulting logLiR a= i 2.5 was adopted 

as a best estimate for the IR luminosity. Table [3] lists the 
estimated total IR luminosities, Lir, with the corresponding 
1 a errors, with and without the uncertainties due to ran- 
dom photometric redshift errors. The adopted approach to 
estimate Lir from 24 /im fluxes has been shown to produce 
SFRs that are in bette r agreement with t h e SFRs determined 
from SED modeling dFranx et al.l l2008t IWuvts et al l 2008) 
and from dust-corrected Ha line fluxes (Muzzin et al. 2010), 
with respect to the often used local luminosity-dependent ap- 
proaches, which can systematically o v er-estimate SFRs by 
a factor of 4-6 dPapovich et alj 120071: iMurphv et all 120091: 
iMuzzin et alj|2010l) . More importantly, our approach adopted 
to estimate Lir from 24 /im fluxes is further supported by the 
detailed anal ysis of the far- I R SED (from 24 fim to 870 /im) 
performed in IMuzzin et alj ((2010) on two ultra-luminous in- 
frared galaxies (ULIRGs) at z ~ 2. 

We convert the Lir to SFR using the iKennicufil i 19981) re- 
lation, SFR(Li R ) = 0.63 Lir/5.8 x 10 9 L^ , where the factor 
of 0.63 is to convert to a pseudo- [Kroup a| d2001|) IMF. The 
estimated SFRs are listed in Table [3] along with the 1 a er- 
rors, with and without the uncertainties due to random pho- 
tometric redshift uncertainties. As shown in Table [3] the 
random uncertainties on Lir and SFRs are dominated by the 
contribution from random photometric redshift errors. Sys- 
tematic template uncertainty (not included in the errors in 
Table [3) can contribute additionally to the total error bud- 
get, with ±0-45 de x variation from logLi R a= 2 5 to logLi R a= \ 
(IWuvts et alJl2008l) . 

The estimated total IR luminosities of the MIPS detected 
sources range from ~ 5.0 x 10 12 L Q to ~ 4.0 x 10 13 L Q , 
with 80% of them having Lir > 10 13 L Q , typical of Hyper- 
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Luminous Infrared Galaxies (HLIRGs), and the remaining 
being ULIRGs. The fraction of HLIRGs in the most mas- 
sive galaxies at 3.0 < z < 4.0 is larger by a factor ~ 10 with 
respect to the fraction of HLIRGs in the IKrieket all (120081) 
sample of massive (M stal - ss 10 11 M^) galaxies at 2 .0 < z < 2.7 
with spectroscopic redshifts (Muzzin et al. 2010). Whereas 
the sample at 2.0 < z < 2.7 is less massive than our mass- 
selected sample at 3.0 < z < 4.0 by a factor of ~ 3, the large 
difference in the fraction of HLIRGs seems to suggest a large 
evolution in the number density of massive HLIRGs from 
z=3.5toz = 2.3. 

As previously noted, the observed 24 /im emission could 
arise from warm/hot dust and PAH molecules, heated by ei- 
ther dust-enshrouded star formation or highly-obscured AGN. 
If all the IR luminosity is associated with dust-enshrouded 
star formation, the resulting SFRs range between ~ 600 to 
~ 4300 M yr" 1 , with the exclusion of Cl-22857, for which 
only an upper limit was derived. These values are tens to sev- 
eral hundreds of times larger than the SFRs estimated from 
SED modeling. On average, the SFRs estimated from the 
24 /im fluxes are a factor of ~ 200 larger than the SFRs es- 
timated from SED modeling. Moreover, three galaxies (Cl- 
6110, Cl-15182, and Cl-23152) have MlPS-derived SFRs of 
the order of 1.3-2.5x 10 3 M Q yr" 1 , whereas the FUV-to-8 /x 
observed SEDs are consistent with zero star formation. This 
suggests that if the 24 /im flux is from star formation, most 
of it must be completely obscured by dust. The MIPS esti- 
mated SFRs translate to specific star formation rates sSFR ss 
10~ 8,8 - 10~ 78 yr" 1 , which would imply that the most massive 
galaxies at z = 3.5 are extremely actively star-forming sys- 
tems that would double their stellar mass in (0.6-7) x 10 8 yrs 
if the derived SFRs were to be sustained at the current lev- 
els. However, little evolution seems to have been found in the 
number density of the most ma ssive galaxies from z = 3.5 to 
z = 1.6 (Marche sini et alj|2009i) . which would imply a growth 
in stellar mass in the most massive galaxies over this redshift 
range by ^30% (> 3 times smaller than the implied growth 
from z = 3.5 to z = 2.5 from the MlPS-derived SFRs), although 
larger evolution would be allowed once systematic uncertain- 
ties are taken into account. Therefore, either the very large 
star-forming activity indicated by the observed 24 /im emis- 
sion has to be very quickly quenched in the majority of the 
most massive galaxies at 3.0 < z < 4.0, or the MlPS-derived 
SFRs are systematically biased by, e.g., contamination from 
AGN activity. 

Indeed, in the local universe, AGN is thought to be the 
dominant source of radiation responsible for the f ar- and mid- 
IR SEDs of galaxies w ith Lir ~ 10 13 L Q (e.g.. iGenzel et alj 
[19981; iLutz et aljfT998h . If the MIR emission is dominated by 
AGN-heated dust, the large fraction of very massive galaxies 
at 3.0 < z < 4.0 with MIPS detection suggests that AGNs are 
very common (>80%) in the most massive galaxies at these 
redshifts. While the fraction of AGNs in dropout galaxies at 
z > 3 is generally estimated to be low (~3 -7%; lSteidel et alj 
120021: lLaird et aljr2006t iReddv et alj|2006l) . estimates of the 
AGN fraction in massive galaxies at 3.0 z < 4.0 are still 
very uncertain. Common AGN activity in massive galax- 
ies has been found a t lower redshifts, wi t h AGN fraction o f 
~ 30% at z ~ 2.5 dPapovich et al.ll2006t IKrieket al.ll2007h . 
Evidence for an increasing fraction of AGN as a func tion of 
stellar mass as been also shown by iRriek et al] (120071) . with 
AGN fraction that could reach as high as 70% for the most 
massive galaxies at 2.0 < z < 2.7. If the observed 24 /im 



emission represents a signature of AGN-heated dust, then our 
results represent further supporting evidence for higher AGN 
fractions at high-z and in the most massive galaxies. 

With the currently available data it is not possible to 
discriminate between dust-enshrouded starburst or highly- 
obscured AGN as the dominant source responsible in heating 
the dust in our sample of very massive galaxies. Significant 
contributions from both sources to the observed MIPS 24 fim 
fluxes cannot be excluded, and their relative importance will 
certainly vary among our sample. However, if most of the 
massive galaxies in our sample have extreme SFRs, as derived 
from the MIPS data, then it is unlikely that we are witnessing 
short-lived bursts, as the duty cycle of the star formation has 
to be long to account for the observed large fraction of MIPS- 
detected sources. This seems to be in contradiction with the 
need for the extreme MlPS-derived star-formation activity to 
be rapidly (< 10 8 years) quenched to avoid overprediction of 
the high-mass end of the SMF of galaxies at z < 3. Moreover, 
the estimated Lir are typical of HLIRGs, for which AGN is 
generally thought to be a significant (if not dominant) source 
of radiation responsible for the very large IR luminosities. Fi- 
nally, for the targeted redshift range, the MIPS 24 /im band 
probes rest-frame wavelengths around ^5.5 /im, where hot 
dust dominates the MIR emission, and the contribution from 
the AGN as the source of the radiation field heating the dust 
becomes increasingly more likely. For all these reasons, the 
very high MlPS-estimated star formation rates are unlikely, 
and we therefore favor AGN instead of starburst activity as 
the dominant source of the observe MIPS 24 /im emission. 

Whatever the source of radiation responsible for heating the 
dust is (dust-enshrouded star formation and/or AGN), the very 
large IR luminosities estimated in our sample of very mas- 
sive galaxies at 3.0 < z < 4.0 show that, despite the already 
very large stellar masses, there is still plenty of gas and dust 
around either the supermassive black holes or the star forming 
regions. 

3.2.1. X- ray em ission 

AGNs can be efficiently identified by their X-ray emis- 
sion, which is thought to be due to up-scattered UV photons 
from the accretion disk. AGN-induced X-ray emission can 
be distinguished from that induced by star formation by the 
hardness ratio and (particularly) the luminosity. Chandra X- 
ray data are available over both the COSMOS and AEGIS 
fields. We have u sed the publicly avai lable X-ray catalogs 
(lLaird et alj 120091 and lElvis et al.ll2009l for the AEGIS and 
COSMOS fields, respectively) to search for X-ray detections 
within our mass-selected sample at 3.0 < z < 4.0. 

Three sources, namely Cl-15182, Cl-19536, and A2- 
15753, are detected in the Chandra images, with hard band 
(2-7 keV), fluxes of (3.8±1.0)x 10" 15 , (1.3±0.2)x 10" 14 , 
and (7.6±1.2)x 10~ 15 erg s" 1 cm" 2 , respe ctively. Assuming 
a po wer-law photon index T = 1.9 (Nand ra & Pounds] 
1994), these fluxes correspond to X-ray luminosi- 
ties L 2 -7kev of (3.8±1.3)xl0 44 , (1.0±0.2)xl0 45 , and 
(5.7±1.5)x 10 44 erg s" 1 , respectively, typical of high- 
luminosity AGNs (L 2 -7keV > 3 x 10 43 erg s _1 )13 Their 
hardness ratios are 0.17±0.3, -0.24±0.11, and 0.31 ±0.1 3 
typica l of narrow-line and high-z obscure d AGN s dBrusa et alj 
120091) . Using Figure 5 in ifreister et ail (120091) . which plots 

18 The quoted errors of L2-7 kcV includes the error due to photometric red- 
shift random uncertainties. 
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the neutral hydrogen column density (A^h) as a function of 
hardness ratio for 15 high-redshift X-ray sources, we esti- 
mated Nh =(4^)x 10 23 , (l!^)x 10 23 , and (6±2)x 10 23 cm" 2 , 
characteristic of moderately obscured AGNs. 

We therefore conclude that the X-ray data for the three de- 
tected sources support the presence of powerful AGNs in all 
three sources, providing further evidence for AGN radiation 
being the dominant source for heating the dust and for the 
large MIPS fluxes. Note however that only high-luminosity 
AGNs with L 2 -7 kev > 10 43 7-43 9 erg s" 1 could have been de- 
tected for sources at z > 3 given the limiting source detec- 
tion depths of the X-ray data over AEGIS and COSMOS. In 
other words, the lack of X-ray detection does not provide in- 
formation on the presence, or lack thereof, of high-luminosity 
AGNs in the remaining 1 1 sources of our mass-selected sam- 
ple. 

3.2.2. Continuum emission from the AGN 

In the previous section, we found that the AGN is likely the 
dominant source of radiation responsible for heating the dust 
and for the large IR luminosities. Moreover, three galaxies 
have X-ray detections with hardness ratios and estimated X- 
ray luminosities typical of obscured high-luminosity AGNs. 
Therefore, the AGN emission could potentially contribute to 
the observed SED, biasing the derived stellar masses. We in- 
vestigate this by subtracting the AGN contribution from the 
observed SED and by re-fitting the corrected SED. Specifi- 
cally, we assume a power-law SED for the AGN, with F v oc 
v a . The value of a has been derived by fitting the rest- frame 
UV and the MIPS 24 /im photometry, with values of a found 
in the range -2.9 < a < -1.7. The maximum AGN contri- 
bution is then set by the rest-frame UV fluxes in combination 
with the 24 /im band, and subtracted from the observed SED. 
The resulting SEDs are finally re-modeled using FAST to de- 
rive stellar masses. We find that the derived stellar masses are 
smaller by typically ^0.08 dex. For two of the three galax- 
ies with X-ray detection, this analysis results in stellar masses 
smaller by only 0.05-0.08 dex, slightly larger than their ran- 
dom errors on M star . For the third galaxy with X-ray detection 
(A2-15753), the AGN contribution is 0.18 dex, the largest in 
our mass-selected sample (although still much smaller than 
the systematic uncertainties due to different SED-modeling 
assumptions). 

We stress that the estimated systematic effects caused by the 
AGN contributions should be strictly regarded as upper lim- 
its, as our approach maximizes, by construction, the contribu- 
tion of the AGN to the observed SED. We therefore conclude 
that these systematic effects are in general small, and cer- 
tainly much smaller than the systematic uncertainties caused 
by the different SED-modeling assumptions and/or by po- 
tential systematic errors in the photometric redshift estimates 
(see § [3). We note that additional contamination of some of 
the medium- and broad-band filter fluxes could be potentially 
caused by the presence of strong AGN line emission. 

4. THE STELLAR MASS FUNCTION AND DENSITY 

4.1. Methodology 

We used the mass-selected sample defined in §|2]to derive 
more accurate measurements of the high-mass end of the SMF 
of galaxies at 3.0 < z < 4.0. T o estimate the observed S MF we 
have followed the analysis in Marchesi ni et alJ J2009), which 
we refer to for a detailed descriptions of the methods used in 
this work. Briefly, we have derived the SMF using two meth- 



ods, the 1 /V max estimator and a parametric maximum likeli- 
hood method. 

For the 1 /Vmax estimator , we u sed the extended version as 
defined bv lAvni & Bahcal] dl980l) . The Poisson erro r in each 
stellar mass bin was computed adopting the recipe of Gehrels 
( 1986). As extensively discussed in the literature, the 1 /V max 
estimator has the advantages of simplicity and no a priori as- 
sumption of a functional form for the stellar mass distribution; 
it also yields a fully normalized solution. However, it can be 
affected by the presence of clustering in the sample. Field-to- 
field variation represents a significant source of uncertainty in 
deep surveys, since they are characterized by small areas and 
hence small probed volumes. The contribution du e to cosmic 
variance to the total error budget is quantified in § 14.21 

We also measured the observed SM F using the STY method 
dSandage. Tammann & Yahill 11979b . which is a parametric 
maximum-likelihood estimator. The STY method has been 
shown to be unbiased with respect to den sity inhomogeneities 
(e.g., [Efstathio u. Ellis & Petersonlll988l). it has well-def ined 
asymptotic error properties (e.g.. lKendall & S tuart 1961]) and 
does not require binning. We have assumed that the number 
density <£>(M) of galaxies is described by a lSchechterl (119761) 
function, 

$(M) = (In 10)^*[10 (M - M * )(1+Q:) ] x exp [- 10^^] , (1) 

where M = log (M star /M©), a is the low mass-end slope, M* = 
log(M* tar /M Q ) is the characteristic stellar mass at which the 
SMF exhibits a rapid c hange in the slope , and 3 >* is the nor- 
malization. Following Marchesin i et al.l (12009b . the best-fit 
solution is obtained by maximizing the likelihood A with re- 
spect to the parameters a and M*. The value of $* is then ob- 
tained by imposing a normalization on the best-fit SMF such 
that the total number of observed galaxies in the sample is 
reproduced. 

4.2. Uncertainties on the Stellar Mass Function 

As discussed in M archesini et al. (2009), uncertainties due 
to small-number statistics, photometric redshift errors, cosmic 
variance, and different SED-modeling assumptions contribute 
to the total error budget of the high-mass end of the high- 
redshift SMF. 

The uncertainties on the SMF due to random photometric 
redshift errors have be en estimated following the recipe in 
Marchesini et aljj20"0 9)PI Specifically, for each galaxy in the 
A'-selected sample, a set of 200 mock SEDs was created by 
perturbing each flux point according to its formal error bar. 
Second, we estimate d the photometric redshift in the same 
way as described in § 12.21 Third, we fit the mock S EDs with 
FAST to estimate stellar masses as described in § 12.31 Fi- 
nally, we have derived SMFs of galaxies with the 1 /V max and 
the maximum likelihood analysis for each of the 200 Monte 
Carlo realizations of the ^-selected sample. The contribution 
to the total error budget of the SMF derived using the 1 /V max 
method due to random photometric redshift errors (cr z ran ) is 
listed in Table [4] and is roughly 0. 1 3 dex , about a factor of 1 .7 
smaller than in Marchesi ni et alJ ( 120091) . Similarly to what 

19 In Marchesini et al. (2009), systematic photometric redshift errors were 
estimated by adopting different template sets to derive the photometric red- 
shifts. Here we decided not to use the other template sets distributed with 
EAZY du e to the significantly w orse resulting Zphai ~ Zspec comparisons, 
whereas in Marchesini et al. 1 2009) they resulted in z p hot — Zspec comparison 
of similar quality, or only slightly worse. 
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found in March esini et alj (12009b . the contribution of random 
photometric redshift errors on the Schechter function param- 
eters of the SMF at 3.0 < z < 4.0 is instead negligible with 
respect to Poisson errors. In fact, Poisson errors largely domi- 
nate the random error budget of the Schechter function param- 
eters due to the complete lack of constraint on the low-mass 
end slope a (see Table[5]l. Because the low-mass end slope is 
not constrained by the NMBS dataset, we have also repeated 
the maximum-likelihood analysis fixing the value of the low- 
mass end slope at a = -1 .0 (corresponding to the value of the 
low-m ass end slope of the S MF of galaxies at 1.3 < z < 2.0 
fromlMarchesin i et alJ2009i) and a = -1 .75 (corresponding to 
the value of the lo w-mass end slope of th e SMF of galaxies at 
2.5 < z < 3.5 from lKaiisawa et al.ll2009h . 

To quantify the uncertainties due to field-to-field varia- 
tions in the dete r minati on of the SMF, we proceeded as in 
iMarchesini et alJ (120071) . Briefly, using the 1 /Vmax method, 
we measured $ ; , the galaxy number density in the stellar mass 
bin AM for the jth field. The contribution to the error budget 
from cosmic variance is estimated using er cv = rm.s(<l> / )/v / 2. 
The final 1 a random error associated with $(M) is then 
a = (cpoi + cr^ + o^nu,) 1 ' 2 , where <7p i is the Poisson error in 
each mass bin. The values of erp i and cr cv are also listed in 
Table H] The contribution to the error budget from cosmic 
variance can also be estimated for a given population using 
predictions from cold dark matter theory and the galaxy bias. 
We have derived the cosmic va riance follow i ng the cosmic 
variance cookbook presented bv lMoster et alj (120101) and us- 
ing the parameters for our survey and for a population of mas- 
sive galaxies withM star > 10 11 M Q , resulting in an uncertainty 
due to cosmic variance of 0.18 dex, in very good agreement 
with our empirical est imate. 

Whereas we refer to M archesini et al.1 (120091) for a complete 
analysis and discussion of the systematic uncertainties due 
to different SED-modeling assumptions, we have repeated 
the whole analysis adopti ng the st ellar pop ulation synthesis 
models of Maraston (120051) with a iKroupal (120011) IMF and 
exponentially declining SFHs with values of the e-folding 
timescale ranging from 100 Myr to 10 Gyr. The resulting sys- 
tematic uncertainties on the SMF measured with the 1 /V max 
method are listed in Table [4] and the corresponding values of 
the Schechter function parameters measured with the maxi- 
mum likeli hood analys is are listed in Table [5] The adoption 
of the Maraston (2005) models result, in general, in smaller 
d erived stellar masses by ~ 0. 15 dex. As previo usly shown 
in lMarchesini et al.1 d2009) and lMuzzin et al.1 (120091) . different 
combinations of adopted metallicity and extinction curve also 
result in systematic effects on the derived stellar masses, al- 
though to a much smaller extent with respect to the biases in- 
troduced by different stellar population synthesis models and 
the specific choices of the adopted SFHs. In particular, dif- 
ferent assumptions on the SFH with respect to those adopted 
in our work (e.g., two-component models of the SFH, or 
exponentially-increasing SFH) can introduce additional sys- 
tematic biases toward both larger and smaller stellar ma sses 
dWuvts et al.ll2007l: iLee et al] |2009: Mar aston et aHHoIol) . In 
§ [5] we consider the systematic uncertainties due to the in- 
clusion of an additional template in the template set used to 
derive photometric redshifts with EAZY. This additional tem- 
plate consists of an old (1 Gyr) and dusty (Ay = 3 mag) single 
stellar population. 

4.3. Stellar Mass Function 
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FIG. 8.— SMFs of galaxies at redshift 3.0 < z < 4.0 from the NMBS (red 
orange, and yellow symbols) and the analysis of Marchesini et al. ( 2009) 
(black and gray symbols). The filled symbols represent the SMF derived with 
the 1/Vmax method, with error bars showing the total 1 a random errors, in- 
cluding photometric redshift errors and field-to-field variations; the red boxes 
include also the systematic unc ertainties due to the different SED-modeling 
assumptions adopted (see j) |2.3> . The solid curves represent the SMFs derived 
with the maximum likelihood analysis, with shaded regions representing the 
1 cr uncertainties. The black error bars and gray shaded area include the sys- 
tematic uncertainties due to different template sets in the photometric redshift 
estimate. The orange hatched area includes also the systematic uncertainties 
due to the different SED-modeling assumptions adopted in our analysis. The 
dotted and dashed b lack curves represen t the predicted SMFs from the semi- 
analytic model of Somerville et al. 12008), where the dashed curve is derived 
from the dotted curve after convolution with a normal distribution of stan- 
dard deviation of 0.25 dex. The NMBS allows us to derive more accurate 
measurements of the high-mass end of the SMF of galaxies at 3.0 < z < 4.0. 



Figure[S]shows the SMF of galaxies at redshift 3.0 < z < 4.0 
derived in this work (colored symbol s) compared to the SMF 
of galaxies at 3.0 < z < 4.0 derived in M archesini et al.l (120091) 
(black and gray symbols). Points with error bars show the 
SMFs derived using the 1 /V max method. The solid curves 
show the SMFs derived with the maximum likelihood anal- 
ysis, while the shaded regions represent their 1 a uncertain- 
ties. The plotted uncerta inties of the SMF measurements from 
Marchesini et ail (120091) . the thick red errors bars, and the yel- 
low shaded area represent the total 1 a random errors, includ- 
ing cosmi c va riance and photometric redshift errors as quan- 
tified in § 14.21 The thin red error bars and the orange shaded 
area include also the systematic uncertainties due to the dif- 
ferent SED-modeling assumptions adopted in this work. 

The large surveyed area (i.e., effec tive area of 0.44 square 
degrees, a factor of ~3 larger than in March esini et al.ll2.Q09t) 
and the accurate photometric redshift estimates allow for the 
determination of the number density of the most massive 
galaxies at 3.0 < z < 4.0 with unprecedented accuracy, as 
clearly shown by Figure [8] f rom the compar i son w ith the 
SMF previously derived by Mar chesini et al.l ((2009). Fig- 
ure [9] shows the comparison between the SMF derived in this 
work and previous measurements of the SMFs of galaxies at 
z ~ 3.5. The high-mass end of the SMF measured in our anal- 
ysis is in good agreement with prev ious measurement s. 

Combined with the results from IMarchesini et al.l d2009l) . 
the number density of the most massive galaxies appears to 
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FIG. 9. — Top panel: Comparison between the SMF at 3.0 < z < 4.0 (yel- 
low and red) derived from this work and previous measurements fr om the lit- 
erature, including the measurement from Marchesini et al. 1 2009 ) (gray and 
black). For the SMFs derived from this work and Marchesini et al. (2009), 
the filled circles represent the measurements using the 1 /V m ax method, while 
the solid curves represent the measurements using the maximum-likelihood 
analysis; the 1 a error bars of the 1 /V max measurements include Poisson 
errors, field-to-field variations, and uncertainties from photometric redshift 
uncertainties (both random and systematic). Similarly for the 1 a error of 
the maximum-likelihood measurements (yellow and gray regions). Previ- 
ous works are plotted as filled stars and d ashed curves lFonta naetalj|200"6t 
F06); open circles and long-dashed curves I Perez-Gonzalez et al. 2008; P08); 
open stars and dot-dashed curves (Eisner et al. 2008; E08); open triangles 
and dotted curves IDrorv et al. 2005; D05). Th e hatc hed green area shows 
the SMF of B-dropout galaxies from Stark et al. (2009). Bottom panel: Sym- 
bols as in the top panel, but now the differences bet ween the SMFs measure - 
ments showed in the top panel and the SMF from Marchesini et al. (2009), 
A<J> = log <P — log<E>£)M09, ale plotted as a function of stellar mass to highlight 
the differences. 
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FIG. 10— R atio of the high-z SMFs ($) and the local SMF ($-m).i; 
ICole et al]|2001l) plotted as function of the stellar mass as measured from 
the maximum-likeliho od analysis. The SMFs a t z ~ 1.6 (blue) and z ~ 2.5 
(green) are taken from Marchesini et al. ( 2009). The shaded regions repre- 
sent the total 1 <r random uncertainties. The vertical dashed and dotted lines 
represent the value of 3 X 10 1 1 Mq and the z = 0. 1 characteristic stellar mass, 



10 Mq, respectively. 



TABLE 5 

Best-fit Schechter function parameters of the SMF 
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NOTE. — The quoted error correspond to the 1 a error estim ated 
from the maximum-likelihood analysis as described in § 14.21 The 
values in parenthesis are the best-fit Schechter function pa rameters 
derived with different SED-modeling assumptions, i.e., the MarastoiJ 
12005) stellar population s ynthes is models, [Kroupa 12001) IMF, solar 
metallicity, [Calzetti et al. 12 000) extin ction curve, and exponentially 
declining SFHs (see § 12.31 and § 14. 2) . Note that the low-mass end 
slope is completely unconstrained. The second and third rows list the 
best-fit Schechter function parameters obtained with fixed a. 



TABLE 4 

SMF AT 3 < Z < 4.0 DERIVED WITH THE 1 /Vmax METHOD 



logM star 
(Mq) 


log<P 
(Mpc~ 3 dex" 1 ) 
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0"Poi 


CTz.ran 


CTcv 


°"MA05,sys 


11.63 
11.48 


-5.282 
-4.842 


+0.390 
-(1.4(19 
+11.2411 
-0.244 


+0.253 
-(1.282 
+0.146 
-0.153 


0.151 
0.117 


0.256 
0.150 


+0.0 
-0.6 
+0 (1 

-0.6 



NOTE. — o = (cP- pm + a\ v + 0'J ra „) 1 / 2 is the total 1 o random error, includ- 
ing the Poisson errors (crp j), the errors due to random photometric re dshif t 
uncertainties (a>, ra „), and the error due to cosmic variance (cr c ,,; see § 14.21 ; 
CMA05.sys is the systematic uncertainty due to the different SED-modeling 
assumptions, i .e., the Maraston (2005) stel lar popu l ation synthesis models, 
IKroupal (2001) IMF, solar metallicity, C alzet ti et a l. 12 000) extinction curve, 
and exponentially declining SFHs (see § |2.3| and § |4.2> . 



have evolved by a factor of ~ 2 from z = 3.5 to z = 2.5, and 
by a factor of ~ 3 from z = 3.5 to z = 1.6. This is shown in 
Figure [10] where the r atio of the high-z SMFs and the local 
SMF from ICole et al.l (1200 ll) is plotted as a function of the 
stellar mass. However, due to the steepness of the high-mass 
end, the implied evolution of the number density translates 
to small growth in stellar mass of the most massive galax- 
ies, by 30%-40% from z = 3.5 to z = 1.6, and by -40% from 
Z = 1 .6 to z = 0. 1, although s ystematic uncertaint i es all ow for 
a larger evolution. Recently, Ivan Dokkum et al] (120101) have 
estimated a growth of a factor of — 2 in the stellar mass of 
massive galaxies from z = 2 to z = 0.1, in apparent contradic- 
tion with our results. However, the selection of the sample of 
massive galaxies in lvan Dokkum etaildlOiOh was very differ- 
ent than ours, as galaxies were selected at a constant number 
density of n = 2 x 10~ 4 Mpc" 3 over the targeted redshift range. 
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This approach selects galaxies with stellar masses logM stal - = 
1 1.45 ± 0.15 at z = 0.1, and logM star = 11.15 ± 0.15 at z = 2, 
a factor of ~ 1.2 and ^2.1 smaller than the typical galaxy 
in our mass-selected sample at z = 3.5. As a consequence 
of the mass-dependent evolution derived in Marchesi ni et alj 
(2009), a smaller growth in the stellar mass of massive galax- 
ies would therefore be derived if the selection were done at 
a value of the number density typical of galaxies with stel- 
lar masses logM star = 1 1 .5 at z = 0. 1 (the average mass of our 
mass-selected sample at z = 3.5). We also note that systematic 
uncertainties in the SED modeling, as well as the choice of the 
z ~ benchmark, can play an important role in the estimate of 
the evolution of the stellar content i n massive galaxies. There- 
fore, we conclude that the results in Ivan Dokkum e tail d20loh 
are broadly consistent with ours once the different selection 
and systematic uncertainties are taken into account. 

Figure [8] also sho ws the SMF predicted from the semi- 
analytic model of Somervil le et alj Q2008) (dotted black 
c urve H This model, built on th e pr evious models describe d 
in lSomerville & Primackl (119991) and lSomerville et al.1 (1200 ll) , 
presents several improvements, including, but not limited to, 
tracking of a diffuse stellar halo component built up of tidally 
destroyed satellites and stars scattered in mergers, galaxy- 
scale AGN-driven winds, fueling of black holes with hot gas 
via Bondi a ccretion, and heating by radio jets. The predic- 
tion from the lSomerville et aT] (120081) semi-analytic model are 
taken from their fiducial WMAP-3 model, which adopts a 
fraction / sca tter = 0.4 of the stars in merged satellite galaxies 
added to a diffuse component distributed in a very extended 
halo or envelope. The dashed curve in Figure [8]represents the 
model-predicted SMF convolved with a normal distribution of 
standard deviation 0.25 dex, intended to rep resent measure- 
ment errors in logM star dFontanot et alJl2009H 21 L 

The comparison between the model-predicted and the 
NMBS-derived SMF provides further supporting evidence for 
the deficit of very massive galaxies at 3.0 < z < 4.0 in the the- 
oretical models of galaxy formation, a disagreement that was 
only marginal with the previously derived SMFs. Without 
the inclusion of the systematic uncertainties, the disagreement 
between the observed and the (convolved) model-predicted 
high-mass end of the SMF of galaxies at 3.0 < z < 4.0 is sig- 
nificant at the 3 a level. The significance of the disagreement 
is reduced to only 1 a if we include the systematic uncertain- 
ties due t o different SED-mod eling assumptio ns as estimated 
in § 14.21 (i.e., adopting t he iMarastoril d2005l) instead of the 
iBruzual & Charlotl (T2003) models, and different SFHs). We 
note that systematic uncertainties due to an evolving IMF can 
play an additional role in reducing the disagreement between 
observed and model-predicted SMFs. 

4.4. Number and Stellar Mass Densities 

The number density, tj, and stellar mass density, p, in mas- 
sive galaxies at 3.0 < z < 4.0 has been estimated by inte- 
grating the SMF at M star > 10 1140 M Q These values are 

20 In the compar ison with the mod el predictions, we decided to use 
only the model of Somerville et al. 12008) as it was the model showing the 
smallest disagreements at the high-mass end in the comparison presented in 
March esini et alj <2009D 

21 We note that the typical random error on the stellar masses for the galax- 
ies in our 3.0 < z < 4.0 sample is smaller than 0.25 dex by a factor of ~ 2, 
due to the combination of accurate photometric redshift estimates and well- 
sampled SEDs delivered by the NMBS 

22 The SMF has been integrated using M sta r = 10' 3 Mq as the upper limit 
of the integral. Due to the exponential behavior of the SMF at the high-mass 



TABLE 6 

Number and stellar mass densities at 3 < z < 4.0 



M star > 10 1L40 M 



M stm > 10 8 M 



log(r? [Mpc- 3 ]) (-6.04^ 2 g) > -5.47 (> -6.00) 

log(p [M Q Mpc- 3 ]) 6.00™;™ (5.49;g 3 °) > 6.09 (> 5.59) 

log (r, [Mpc- 3 ]) -5.55™|5 (-6.04+019) -3.61+°^ (-3.79!|- 33 ) 

log(p[M Mpc- 3 ]) 6.00™;! 2 (5.48;0 ; |5) 6.74;° ; 3 > (6.45™« 2 ) 



log (r, [Mpc- 3 ]) -5.55™ ] 8 (-6.05 + 0- 2 0) -1.90+°'' (-2.06!' 28 ) 

log(p[M Mpc- 3 ]) 6.01^ |4 (5.48+g Jg) 7.25+» 3| (7.04^ ™) 



NOTE. — Number density, rj, and stellar mass density, p, at 3.0 < z < 
4.0 estimated by integrating the best-fit Schechter SMF over the spec- 
ified stellar mass range. The quoted 1 a errors include Poisson errors, 
errors due to photometric redshift uncertainties, and errors due to cos- 
mic variance. The values in parenthesis are the resu lts corresponding 
to the different SED-modeling assu mptions, i.e., the Marastoq <2005H 
stel lar population synt hesis models, Kroupa (2001) IMF, solar metallic- 
ity, Calzetti etal. (20 00) e xtinction curve, and exponentially declining 
SFHs (see § !2.3l and § 14.21 . The third and fifth rows list the number den- 
sities estimated fixing the low-mass end slope at a = —1 .0 and a = — 1 .75, 
respectively; the fourth and sixth rows list the stellar mass densities esti- 
mated fixing the low-mass end slope at a = —1.0 and a = —1.75, respec- 
tively. 

listed in Table [6] together with the corresponding total 1 a 
errors. Table [6] also lists the 1 a lower limits of the num- 
ber and stellar mass densities of galaxies more massive than 
10 8 Mm, as well as the densities estimated with the second 



set of SED-modeling assumptions adopted (see § 12. 3t . The 
number and stellar mass densities estimated by fixing the low- 
mass end slope at a = -1.0 and a = -1.75 are also listed 
in Table [6] Compared to the total stellar mass density in 
galaxies with M. t a r > 10 10 M Q at 3.0 < z < 4.0 estimated by 
Marchesini et al] (120091) . the contribution of the galaxies in 
our mass selected-sample is ~ 8^ 3 % (not including system- 
atic errors). Note however that this estimate depends strongly 
on the value of the l ow-m ass end slope of the SMF derived 
in iMarchesini et al.l (120091). which i s still very poorly con- 



strained (e.g., Kajisawa et al. 2009; March esini et 



orlycon- 
aiT l2009l: 



iReddv & Steidelll2009l) 

In recent years, there have been several claims of the ex- 
istence of a population of massive (and evo lved) galaxies at 
even larger redshifts, i.e., z > 4 — 5 (e .g., lYan et al] 120061 : 
IWiklindet al.ll2008UMancini et alj|2009l) . 

In particular, Wiklind et al. (2008) reported of a significant 
population of massive galaxies at 4.9 < z < 6.5 found in the 
125 arcmin 2 GOODS-South field dominated by old stellar 
populations with M stal = (0.3-3) x 10 11 M . In their sam- 
ple there are only two o bject with M sfar > 2.5 x 10 11 M© 
(one already identified by Mobash er et al.ll2005l as a candi- 
date for a massive, evolved galaxy at z ~ 6.5), implying a 
stellar mass density at z ~ 5.7 of p(M sVdr > 10 UA0 M Q ) « 
2 x 10 6 M© Mpc" 3 , and suggesting no evolution of the stellar 
mass density in the most massive galaxies over t he ~ 800 Myr 
interva l from z ~ 5.7 to z = 3.5. However, iDunlop et al.1 
(120071) have concluded that there is no convincing evidence 
for any galaxy with stellar mass M stal - > 2 x 10 11 M and 
z > 4 in the GOODS-South field, which implies a much 

end, the estimated stellar mass density does not depend significantly on the 
specific value of this limit. 
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stronger evolution of the stellar mass density i n very massive 
galaxi es in the first 1.5 Gyr of the universe. Wikli ndet alj 
(2008) also estimated a stellar mass density in galaxies with 
M st3 _ T > 10 10 - 8 M Q of p > 5 x 10 6 M Mpc" 3 (after correction 
for the different IMF). Combini ng the results from our analy- 
sis and Marchesini et al. (2009), it implies an increase of the 
stellar mass density in massive galaxies by a factor of ~ 1.6 
fromz ~ 5.7 to z = 3.5. 

SMFs and stellar mass densities have also been estimated 
for V- and i-dropout galaxies at z ~ 5, and z ~ 6, respec- 
tively dYan et all 120061: lEvles et all 120071; IStark et all 120071; 
IMcLure et all I2009t IStark et al.l 120091) . A comparison with 
these results is not straightforward. First, these are all 
optically-selected samples, which can be potentially biased 
against massive and evolved galaxies. In contrast, our sam- 
ple is a mass-complete sample constructed from a ^-selected 
catalog. Second, the stellar mass ranges probed by these stud- 
ies are very different from ours. Our mass-selected sample 
probes the most massive galaxies, i.e., those with M star > 
2.5 x 10 11 Mq, whereas all the above works probe galax- 
ies with typically much smaller stellar masses, i.e., in the 
range 1 9 — 1 1 1 M Q . A rough estimate of the evolution 
of the stellar mass density in galaxies more massive than 
2 x 10 9 M can be performed by comparing the stellar mass 
densities derived by the above studies with the stellar mass 
density obtained by co mbining our results with those from 
iMarchesini et alj (120091) at lower stellar masses, and extrapo- 
lating the Schechter function to the stellar mass limit probed 
by the above studies. After correcting for the different IMFs, 
the stellar mass density in galaxies more massive than 2 x 
10 9 Mm atz~5 and z~ 6 is p= ■ (2.3-6.3) x 10 6 M^ Mpc" 3 
( IMcLure et al.1120091; IStark et al.ll2007l; IStark et al.ll2009l). and 
p= (Q .7-4.3) x 10 6 M W Mpc" 3 dYan et al.l|2006UEvleset all 
120071; IMcLure et al.ll2009t IStark et al.ll2Q09h . respectively. At 
z = 3.5, we estimate p = 1.5 x 10 7 M Q Mpc" 3 , which im- 
plies an evolution of the stellar mass density in galaxies with 
M stal - > 2 x 10 9 M by a factor of - 2 - 7 and - 3 - 22 from 
z ~ 5 and z ~ 6 to z = 3.5, respectively. We stress however 
that the estimated evolution from z > 4 is very uncertain, and 
affected by both large uncertainties in the SMF of galaxies at 
z > 4 (especially at the high-mass end), as well as poor con- 
straints on the low-mass end slope of the SMF of galaxies at 
3.0<z<4.0. 

5. MASSIVE, OLD AND DUSTY GALAXIES AT 2 < Z < 3? 

In this section we consider the case of adding an additional 
template to the EAZY template set used to derive photometric 
redshifts. 

Previous searches for old and massive galaxies at z > 4 
highlighted the difficulty in unambiguously identifying old 
and massive objects at extrem e redshifts on t he basis of spec- 
tral fitting. In particular. iDunlop et al.l ([2007) have shown that 
equally acceptable solutions could be obtained at z ~ 5 with 
high stellar masses (M stal ^6xlO n M Q ) and low extinction 
(Ay ~ 0.4 mag) and at z ~ 2 with moderate stellar masses 
(Af stal ~7x 10 10 M Q ) and high extinction (A v ~ 3.8 mag). 

To robustly estimate photometric redshifts, the template set 
needs to be large enough that it spans the broad range of multi- 
band galaxy colors and small enough that the color and red - 
shift degeneracies are kept to a minimum (e.g.,|Benitej2000). 
The default template set used in this work was carefully 
constr ucted and tested in iBrammer. van Dokkum. & Coppil 
(2008). It has been shown to satisfy the requirements for a 
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FIG. 1 1 . — SMFs of galaxies at 3 < z < 4.0 measured by adopting the ad- 
ditional old and dusty template in the estimate of the photometric redshifts. 
The gray hatched region, and the black solid line and empty circles represent 
the SMF measured adopting the default EAZY template set. The colored re- 
gions represent the SMFs derived by adopting the additional old and dusty 
template in the estima te of the photometric reds hifts. The green region is 
obtained by using the Bruzual & Chariot (2003) stellar population synthe- 
sis mo dels, while the hatched red region is obtained by using the IMarastorl 
12005) models. The dotted and dashe d curves represent the predicted SMFs 
from the semi-analytic model of Somerville et al. 1 2008). 

satisfactory template set, providing significantly reduced sys- 
tematic effects and smaller scatter in the z p h t versus z spe c at 
all redshifts. This template set already includes a dusty star- 
burst model (50 Myr old and Ay = 2.75 mag). Here we in- 
clude an additional template representative of an old (1 Gyr; 
t = 100 Myr) and very dusty ( Ay = 3 mag) galaxy, simi - 
lar to the reddest template used in Blan ton & Roweisl d2007l) . 
and we repeat the whole analysis. We note that in the local 
Universe old stellar populations are usually not (very) dust- 
obscured, and that it remains to be seen whether this template 
is physically plausible for the high-redshift Universe. 

For eight galaxies (CI -4890, CI -6 110, CI -7340, Cl- 
15367, Cl-18825, A2-6835, A2-15753, and A2-18070), the 
resulting photometric redshift estimates formally lie at z < 3. 
The three objects in the AEGIS field have now z ~ 2.9, hence 
just below our redshift selection window, consistent with their 
redshift probability functions and with A2-15753 having rest- 
frame UV color typical of an LBGs. The five objects in the 
COSMOS field are instead shifted to much lower redshifts, 
i.e., z ~ 2.4, significantly lower than allowed for by their red- 
shift probability functions plotted in FigureQ] 

We used FAST to refit stellar population synthesis models 
to the eight galaxies that moved to z < 3. Six of the eight 
are best fitted by an old stellar population (age«2 Gyr, i.e., as 
old as the age of the universe at their redshifts), large stellar 
masses (M star « 2xlO n M ), and large values of extinction 
(Ay « 2.1 mag). The remaining two objects are instead best 
fitted by very dusty (Ay « 3.2 mag), young (age«50 Myr) 
starbursts. 

The SMFs of galaxies at 3.0 < z < 4.0 measured by adopt- 
ing the additional old and dusty template in the estimate of 
the photometric redshifts are shown in Figure [TT| As shown 
in this figure, the systematic effect on the derived SMF of 
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galaxies at 3.0 < z < 4.0 due to the inclusion of this addi- 
tional (although generally not included) template is signifi- 
cant, larger than the systematic effect caused by the stellar 
population synthesis model assumption alone, and bringing 
the observed SMF of galaxies at 3.0 < z < 4.0 in better agree- 
ment with the SMF predicted from semi-analytic models. 

We note however that the EAZY and FAST best-fit models 
using the additional old and dusty template do not provide sta- 
tistically better modeling of the observed SEDs with respect 
to the default template set without the old and dusty template. 
This further shows the ambiguity and difficulty in character- 
izing the population of massive galaxies at z > 3, even with 
our adopted dataset, which is characterized by exquisite wave- 
length coverage from the ultra-violet to the mid-infrared. 

In summary, our sample of massive galaxies at 3 < z < 4 
could potentially be contaminated (up to ^50%) by a pre- 
viously unrecognized population of massive, old, and very 
dusty galaxies at z ~ 2.6. We note that the existence of such a 
population would be an important and puzzling result in itself. 

6. SUMMARY AND CONCLUSIONS 

In this paper we have used the far-ultraviolet to mid-infrared 
coverage of the NMBS to derive the observed and rest-frame 
properties of a complete sample of galaxies at 3.0 < z < 4.0 
with M s ta r > 2.5 x 10 11 Mq, and to provide more accurate 
measurements of the high-mass end of the SMF of galax- 
ies at 3.0 < z < 4.0. With the addition of five medium- 
bandwidth NIR filters, NMBS delivers accurate photometric 
redshift, cr z /(l +z) ~ 0.02, for Zf-selected sources at z > 1.5, 
and provides well-sampled SEDs in the critical wavelength 
regime around the Balmer/4000 A breaks, allowing, for 
the first time, the accurate detection of the Lyman and the 
Balmer/4000 A breaks simultaneously. Combined with its 
large surveyed area, ~ 0.5 squared degree, it allowed us to 
construct a statistical significant and representative sample of 
14 very massive (M stal - > 2.5 x 10 11 M ) galaxies over the 
redshift range 3.0 < z < 4.0. 

The typical very massive galaxy at 3.0 < z < 4.0 is red and 
faint in the observer's optical, with a median r-band magni- 
tude of (r tot ) = 26.1. The median H-K color is 1.2, with only 
two object having H-K < 0.9 at a significant level, highlight- 
ing the efficiency of the H-K color technique in selecting 
gala xies at z > 3 with prominent bre aks in the rest-frame opti- 
cal dBrammer & van Dokk um 2007). The median rest-frame 
U — V color (U -V) = 1.6 is similar to local Sb spiral galax- 
ies, although we find a range in U — V colors, from the typical 
color of nearby irregular to those of local elliptical galaxies. 
The median U — V slope is (f3) = -0.36, indicating a relatively 
flat spectrum in F\. Intriguingly enough, the distribution of 
UV slopes of the mass-selected sample at 3.0 < z < 4.0 is very 
different from the distributions of UV slopes of UV-selected 
galaxies at z > 2.5 as well as the H-K galaxies at z ~ 3.7 
so far discovered, which show distributions of j3 peaked at 
f3 < -1 .6. This difference is most likely due to the very differ- 
ent ranges in stellar mass prob ed by the different samples with 
the H-K > 0.9 galaxies in lBrammer & van Dokkuml (120071) 
being on average a factor of ^15 less massive than our sam- 
ple, and the typical UV-sel ected galaxies having masses in the 
range 10 9 - 10 11 M (e.g. Shap ley etalJ l2001: MagdiTetaD 
l2010h . 

By constructing a mass-limited sample from a /f-selected 
catalog with accurate photometric redshifts rather than the 
typical color-selection techniques, we were able to find a pop- 



ulation of galaxies mostly complementary to the typical pop- 
ulation of dropout galaxies at z ~ 3-4. Specifically, we have 
shown that only 57% of the mass-selected sample have ob- 
served optical colors that satisfy either the U- or the B-dropout 
color criteria. However, ^50% of these galaxies are too faint 
in the observed optical to be included in typical spectroscopic 
samples of LBGs. 

From the SED modeling, our complete sample of massive 
galaxies at 3.0 < z < 4.0 seems to show a range in stellar pop- 
ulation properties. About 40%-60% of the sample is charac- 
terized by ages consistent with the age of the universe at the 
targeted redshifts, suggesting that the bulk of the stellar mass 
in these systems was formed at very early times. Dust seems 
to be quite ubiquitous in massive galaxies at 3.0 < z < 4.0, 
with a median extinction of (Ay) = 1.0 mag. About 30% of 
the sample have SFR estimates from SED modeling consis- 
tent with no star formation activity, while the rest of the sam- 
ple is characterized by significant star formation activity, as 
high as several hundreds solar masses per year. Of particular 
interest is the z = 3.54 galaxy CI -22857, which has an esti- 
mated stellar mass of ~ 3 x 10 11 M Q , a maximally-old age, 
and completely suppressed star formation. This galaxy is also 
not detected in the Spitzer-MlPS 24 //m data, further support- 
ing its quiescent nature. Recent ultra-deep NIR spectroscopic 
observations have confirmed a massive galaxy at z = 2. 1 8 with 
strongly suppressed star f ormation, pushing back its forma- 
tion redshift to z > 4-7 dKriek et all 120091) . Spectroscopic 
confirmation of the quiescent nature of Cl-22857, as well as 
other objects in our mass-selected sample, is of paramount 
importance, as it would provide even stronger evidence that 
massive galaxies formed their stars extremely efficiently very 
early in time. 

Surprisingly, most (> 80%) of the massive galaxies 3.0 < 
z < 4.0 are detected in the MIPS 24 /im data. The total 
IR luminosities estimated from the observed 24 /mi fluxes 
range from 5 x 10 12 to 4.0 x 10 13 L , typical of ULIRGs and 
HLIRGs, implying extreme dust-enshrouded star-formation 
rates (^600-4300 M yr" 1 , tens to several hundreds of times 
larger than the SFRs estimated from SED modeling), or very 
common heavily-obscured AGNs, or both in the most mas- 
sive galaxies at z = 3.5. Whereas it is not possible to discrim- 
inated between AGN or starburst as the dominant source re- 
sponsible in heating the dust, we favor AGN as a significant, 
if not a dominant contributor. Specifically, the reasons for this 
are trifold. First, the extreme MlPS-derived SFRs cannot be 
sustained for more than ~ 10 s yrs without overprediction of 
the high-mass end of the SMFs of galaxies at z < 3. This 
seems in contradiction with the very large fraction of MIPS 
detections, which implies long duty cycle of the star forma- 
tion. Second, 80% of the MlPS-detected sources are HLIRGs. 
In the local universe, AGN is thought to be the dominant 
source of radiation responsible for the far- and mid-IR SEDs 
of HLIRGs. Third, for the targeted redshift range, the 24 /im 
band probes the rest-frame wavelengths from ~ 4.8 /im to 
^ 7.1 /im, where hot dust dominates the MIR emission, and 
the contribution from an AGN as the source of the radiation 
field heating the dust becomes increasingly more likely. If 
the MIPS-24 /im emission is dominated by AGN-heated dust, 
the large fraction of very massive galaxies at 3.0 < z < 4.0 
with MIPS detection suggests that AGNs are very common, 
providing further supporting evidence for the coevolution of 
massive galaxies and AGN. We note that three galaxies are 
detected in the X-ray, with 2-7 keV luminosities and hardness 
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ratios typical of obscured, high-luminosity AGNs. Observa- 
tions at longer wavelengths (e.g., in the far-IR), as well as 
other signatures of AGN (e.g., detection of narrow-emission 
lines) are necessary to constrain the occurrence of AGNs in 
this sample and to discriminate between dust-enshrouded star 
formation and heavily-obscured AGN. 

We have significantly improved the measurements of the 
high-mass end of the SMF of galaxies at 3.0 < z < 4.0 (the 
complete analysis of the evolution of the SMF of galaxies over 
the redshift range 0.5 < z < 4.0 from NMBS will be presented 
in Marchesini et al., in prep.). The accurate photometric red- 
shifts and the large surveyed area allowed us to significantly 
reduce the contributions of photometric redshift errors and 
cosmic variance to the total error budget. The measured high- 
mass end is in very good agreement with previous measure- 
ments, providing further supporting evidence for the existence 
of a significant number of very massive galaxies at z = 3.5. 
Combined with the results from Marchesini et al. (2009), the 
number density of the most massive galaxies appears to have 
evolved little from z = 3.5 to z = 1.6, with a larger subsequent 
evolution down to z ~ 0.1. These results are broadly consis- 
tent with the growth of a factor of ~ 2 in the stellar mass of 
massive galaxies from z = 2 to z = 0.1 recently estimated in 
Ivan Dokkum et all (12010). once the different sample selection 
and systematic uncertainties are taken into account. The con- 
tribution of M star > 2.5 x 10 n M© galaxies to the total stellar 
mass budget at 3.0 < z < 4.0 in galaxies with M star > 10 10 M 
is ~ 8±3 3 %, although this estimate strongly depends on the 
value of the low-mass end slope of the SMF, which is still 
very poorly constrained. The stellar mass density in galaxies 
more massive than 2 x 10 9 M Q seems to evolve by a factor 
of 5 ± 3 and 13 ± 10 from z ~ 5 and z ~ 6, respectively, to 
Z = 3.5. 

Our measurement of the high-mass end of 3.0 < z < 4.0 
seems to exacerbate the disagreement between the observed 
number densities of massive galaxies and those predicted 
by the latest generatio n of galaxy formation models (e.g., 
Some rville et al.1 120081) . The disagreement between the ob- 
served and the model-predicted high-mass end of the SMF at 
3.0 < z < 4.0 is significant at the ~3 u level if only random 
errors are considered. However, systematic errors dominate 
now the total error budget at 3.0 < z < 4.0, leading to uncer- 
tainties of a factor of ~ 8 in the densities at the high-mass end. 
When systematic uncertainties due to different SED-modeling 
assumptions are included, the found disagreement between 
observed and model-predicted SMFs is only marginally sig- 
nificant. We finally note that additional systematic uncertain- 
ties on the high-mass end of the 3.0 < z < 4.0 SMF could be 
potentially introduced by either 1) the intense star-formation 
activity and/or the very common AGN activity as inferred 
from the MIPS 24 /im detections, and/or 2) contamination by 
a significant population of massive, old, and dusty galaxies 
at z ~ 2.6 previously unrecognized. This might indicate that 
the high-mass end of the SMF cannot be properly constrained 
without further spectroscopic data. 

The NMBS has allowed us to study with unprecedented ac- 
curacy the population of very massive galaxies at 3 .0 < z < 



4.0, thanks to its wide surveyed area, the accurate photometric 
redshifts, and the well-sampled SEDs in the rest-frame opti- 
cal. To further improve the characterization of the galaxy pop- 
ulation at the high-mass end at 3.0 < z < 4.0, it is necessary to 
significantly increase the sample size, which is currently com- 
prised of only 14 sources. NMBS-II, a shallow-wide accepted 
NOAO Survey Program specifically designed to further con- 
strain the population of very massive galaxies at z > 2, will 
image an area of sky a factor of ~ 10 larger than NMBS, re- 
sulting in a significant increase in the number of very massive 
galaxies out to z ~ 3 .5 with accurate photometric redshifts and 
well-sampled SEDs. 

Follow-up multi-object spectroscopic observations in both 
the optical and in the NIR are of vital importance to confirm 
the redshifts and to better characterize the properties of the 
most massive galaxies at z = 3.5, including AGN and Lyman- 
a emitter fractions, AGN and/or starburst contamination of 
the optical-to-MIR SEDs, superwind outflows, star formation 
rates, and mass-to-light ratios. However, to probe the rest- 
frame wavelength regime red-ward of ~ 5000 A, and hence 
to robustly constrain the star formation histories and SFRs, 
measure metallicities, absorption lines and velocity disper- 
sions from rest-frame optical features, will require NIRSPEC 
on the James Webb Space Telescope. The estimated total IR 
luminosities typical of HLIRGs make the very massive galax- 
ies at 3 < z < 4.0 in our sample ideal candidates for follow- 
up observations with the Atacama Large Millimeter Array 
(ALMA). ALMA will be crucial in constraining the amount 
of dust and gas in these systems, as well as discriminating 
between dust-enshrouded star-formation and obscured AGN 
activity. Moreover, it will also allow for measurements of the 
kinematics in these systems, providing an independent esti- 
mate of the dynamical masses of the most massive galaxies at 
3.0<z<4.0. 

Finally, to fully characterize the population of galaxies at 
3 < z < 4, the analysis performed in this work has to be ex- 
tended to lower stellar masses. This will necessarily require 
very deep imaging with NIR medium band-width filters to 
provide very accurate photometric redshifts and well-sampled 
SEDs down to faint /f-band magnitudes. 
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